#include "../Okada.hxx"

const double Okada::eps(1e-6);

Okada::Okada(const double &lambda, const double &mu, const double &slip,
             const double &Opening,
             const double &Width,
             const double &Length,
             const double &Strike,
             const double &Dip,
             const double &Rake,
             const FTensor::Tensor1<double,3> &Origin):
  alpha((lambda+mu)/(lambda+2*mu)), width(Width),
  length(Length),
  alp1((1-alpha)/2),
  alp2(alpha/2),
  alp3((1-alpha)/alpha),
  alp4(1-alpha),
  alp5(alpha),
  dislocation(slip*cos(radians(Rake)),
              slip*sin(radians(Rake)),Opening),
  origin(Origin),
  rotation_strike(cos(radians(Strike)),sin(radians(Strike)),0,
                  -cos(radians(Strike)),0,-1)
{

  double dip(radians(Dip)), cos_dip(cos(dip)), sin_dip(sin(dip));
  if(std::abs(cos_dip)<eps)
    {
      cos_dip=0;
      if(sin_dip>0)
        sin_dip=1;
      else
        sin_dip=-1;
    }
  cos_dip2=cos_dip*cos_dip;
  sin_dip2=sin_dip*sin_dip;
  sin_cos_dip=sin_dip*cos_dip;
  rotation_dip(0,0)=1;
  rotation_dip(0,1)=rotation_dip(0,2)=rotation_dip(1,0)=rotation_dip(2,0)=0;
  rotation_dip(1,1)=cos_dip;
  rotation_dip(1,2)=-sin_dip;
  rotation_dip(2,1)=sin_dip;
  rotation_dip(2,2)=cos_dip;

  z_transform_dip(0,0)=1;
  z_transform_dip(0,1)=z_transform_dip(0,2)=0;
  z_transform_dip(1,1)=cos_dip;
  z_transform_dip(1,2)=-sin_dip;
  z_transform_dip(2,2)=-cos_dip;
}
